Fusion de imagenes

OBJETIVO PRINCIPAL

La práctica se centra en el Parque Natural de Punta Entinas del Sabinar, una pequeña reseña del lugar la podemos encontrar en el siguiente enlace:

https://www.juntadeandalucia.es/medioambiente/portal/areas-tematicas/espacios-protegidos/legislacion-autonomica-nacional/reservas-naturales/reserva-natural-punta-entinas-sabinar.

El proposito de esta práctica es hacer una fusion con las imagenes de Landsat. Por un lado, las imagenes mutiespectrales y por el otro la pancromática de resolución menor (15 metros).

CARGAMOS LIBRERIAS Y DATOS DE LA ZONA DE ESTUDIO

library(raster)

## Loading required package: sp

library(terra)

## terra 1.7.3

library(rgdal)

## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-4, (SVN revision 1196)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/javiv/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:\Program Files\PostgreSQL\13\share\contrib\postgis-3.0\proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

## 
## Attaching package: 'rgdal'

## The following object is masked from 'package:terra':
## 
##     project

setwd("C:/MASTER_GOFOREST/2023/SENSORES/tarea_isable_2")



LST <- list.files(".", recursive = TRUE, full.names= TRUE, pattern = "B[12345678].TIF")
LST

## [1] "./LC08_L1TP_200034_20230228_20230228_02_RT_B1.TIF"
## [2] "./LC08_L1TP_200034_20230228_20230228_02_RT_B2.TIF"
## [3] "./LC08_L1TP_200034_20230228_20230228_02_RT_B3.TIF"
## [4] "./LC08_L1TP_200034_20230228_20230228_02_RT_B4.TIF"
## [5] "./LC08_L1TP_200034_20230228_20230228_02_RT_B5.TIF"
## [6] "./LC08_L1TP_200034_20230228_20230228_02_RT_B6.TIF"
## [7] "./LC08_L1TP_200034_20230228_20230228_02_RT_B7.TIF"
## [8] "./LC08_L1TP_200034_20230228_20230228_02_RT_B8.TIF"

LST <- lapply(1:length(LST), function (x) {raster(LST[x])})

#cargamos en un stack las bandas multiespectrales y por otro lado 
#la banda Pancromatica

B8<-LST[[8]]
LST<-LST[1:7]
LST_stack <- stack(LST)

# we do a clip centered in our region of interest in this case 
#is the Natural Park "Punta entinas del Sabinar" located in El Ejido, Almería.

PES<-readOGR("./PES.shp")

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\MASTER_GOFOREST\2023\SENSORES\tarea_isable_2\PES.shp", layer: "PES"
## with 1 features
## It has 6 fields

RECORTAMOS EL ÁREA DE ESTUDIO

En este caso se trata del Paque Natural de Punta Entinas del Sabinar, unicado entre los municipios de El Ejido y Roquetas de Mar.

#we could not do a stack of all of our images because Pancromatic band
#has less scale or size of pixel compared with the others, the multiespectral images.

PES_stack_crop <- terra::crop(LST_stack, PES)
plot(PES_stack_crop)

B8<-terra::crop(B8,PES)
plot(B8)

# 1- DATA REESCALE

PES_15<-resample(PES_stack_crop,B8)

dim(PES_15)

## [1] 375 851   7

dim(B8)

## [1] 375 851   1

CALCULO DE ESTADISTICAS : LA MEDIA Y LA DESVIACION TIPICA

Procedemos a calcular las estadisticas media y desviacion tipica de cada una de las bandas multiespectral.

#Calculate mean and standar desviation of each band of multiespectral images

#extract values of each band of image

sts <- data.frame(B_mean = numeric(),B_sd = numeric())

for (i in 1:ncol(PES_15@data@values)){
  sts[i,1] <-mean(PES_15@data@values[,i])
  sts[i,2] <-sd(PES_15@data@values[,i])
  rownames(sts)[i]<-paste0("B",i, sep="")
}

#añadimos las estadisticas mean y sd de la banda pancromatica a tabla de estadisticas

sts[c("B8"),1]<-c(mean(B8[]))
sts[c("B8"),2]<-c(sd(B8[]))

APLICAMOS EL METODO TORUS

#Transformamos la informacion de los píxeles de la imagen numericamente con la Fórmula de Torus.
# Calculate the parameters of the transformation image: Torus's Method 

#we do the same loop to calculate ai and bi from each spectral band

# ai =  sd_Bi/sd_B8_PAN

# bi =  mean_Bi-((sd_Bi/sd_B8_PAN)*mean_B8_PAN)

# And finally:

#PAN_Bi = ai * PAN +bi

# WE CALCULATE

#define a data.frame with a and b like empty comlumns
rm(param)

## Warning in rm(param): objeto 'param' no encontrado

param<-data.frame(a= numeric(),
                  b= numeric())

#calculate each index for each spectral band
for (i in 1:ncol(PES_15@data@values)){
  param[i,1] <-sts$B_mean[i]/sts$B_mean[8]
  param[i,2] <-sts$B_mean[i]-((sts$B_sd[i]/sts$B_sd[8])*sts$B_mean[8])
  rownames(param)[i]<-paste0("B",i, sep="")
}

#PAN_B1 = a1 * B8 +b1

rm(PAN_stk)

## Warning in rm(PAN_stk): objeto 'PAN_stk' no encontrado

PAN_stk<-list()

for (i in 1:ncol(PES_15@data@values)){
  PAN_stk[i] <-param$a[i]*B8 +param$b[i]
  names(PAN_stk[[i]])<-paste0("B",i, sep="")
}

## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = param$a[i] * B8 + param$b[i]): implicit
## list embedding of S4 objects is deprecated

PAN_stk<-stack(PAN_stk)
plot(PAN_stk)

DEFINIMOS EL FILTRO POR LO BAJO

#Define the filter (kernel1)
# Apply the filter that will smooth each image

kernel1 <- matrix(c(0.0039060000, 0.015625000, 0.023438000, 0.015625000, 0.0039060000,
                    0.015625000, 0.062500000, 0.093750000, 0.062500000, 0.015625000,
                    0.023438000, 0.093750000, 0.14062500, 0.093750000, 0.023438000,
                    0.015625000, 0.062500000, 0.093750000, 0.062500000, 0.015625000,
                    0.0039060000, 0.015625000, 0.023438000, 0.015625000, 0.0039060000), 
                  nrow = 5)

print(kernel1)

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 0.003906 0.015625 0.023438 0.015625 0.003906
## [2,] 0.015625 0.062500 0.093750 0.062500 0.015625
## [3,] 0.023438 0.093750 0.140625 0.093750 0.023438
## [4,] 0.015625 0.062500 0.093750 0.062500 0.015625
## [5,] 0.003906 0.015625 0.023438 0.015625 0.003906

#Apply the filter to all multiespectral images transformed
filtered<-list()
for (i in 1:ncol(PES_15@data@values)){
  filtered[i] <-focal(PAN_stk[[i]],kernel1)
  names(filtered[[i]])<-paste0("B",i, sep="")
}

## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = focal(PAN_stk[[i]], kernel1)): implicit
## list embedding of S4 objects is deprecated

filtered_stack<-stack(filtered)
plot(filtered_stack)

APLICAMOS EL FILTRO POR LO ALTO

#When we alpply the low filter o the images next step is subtract this images with originals images


Detail<-list()

for (i in 1:ncol(PES_15@data@values)){
  Detail[i] <-PAN_stk[[i]]-filtered[[i]]
  names(Detail[[i]])<-paste0("DetB",i,sep="")
}

## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PAN_stk[[i]] - filtered[[i]]): implicit
## list embedding of S4 objects is deprecated

Detail_stk<-stack(Detail)
plot(Detail_stk)

### SUMA DE BANDAS: DETALLE + ORIGINAL

#once we applied the "filtro por lo alto" we sum the bands 


fusion<-list()

for (i in 1:ncol(PES_15@data@values)){
  fusion[i] <-PES_15[[i]]+Detail_stk[[i]]
  names(fusion[[i]])<-paste0("FusB",i,sep="")
}

## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated

## Warning in `[<-`(`*tmp*`, i, value = PES_15[[i]] + Detail_stk[[i]]): implicit
## list embedding of S4 objects is deprecated

fusion_stk<-stack(fusion)
plot(fusion_stk)

VISUALIZACION DE LA IMAGEN FUSIONADA FINAL

Representamos las bandas del infrrarrojo, del rojo y el verde.

plotRGB(fusion_stk, r=5, g=4, b=3, scale=maxValue(fusion_stk[[5]]), stretch= "hist")

#writeRaster(fusion_stk,"./stack_fus_automatic.tiff")

COMPARACION DE UNA SEGMENTACION EN LAS IMAGENES FUSIONADAS Y NO FUSIONADAS

#una vez que esta echa la funsion de las imagenes puedo hacer una clasificaion no supervisada y
#comparar si existen variaciones en la clasificacion gracias a la fusion de la imagen

b<-fusion_stk
b[is.na(b)]<-0
kMeansResult <- kmeans(scale(b[]), centers=6)
kMeansResult$centers

##        FusB1      FusB2      FusB3       FusB4      FusB5      FusB6      FusB7
## 1  0.3163838  0.3307411  0.3548348  0.37616767  0.4628357  0.2950999  0.1384995
## 2 -0.6618898 -0.6953083 -0.6977723 -0.62983171 -0.3802437 -0.2323143 -0.1768918
## 3  1.7718918  1.7941766  1.7616567  1.71606633  1.4940061  1.4701645  1.4996942
## 4  1.0441380  1.0686035  1.0752250  1.07383480  1.0003096  0.8885167  0.8020898
## 5 -1.0160851 -1.0339415 -1.0872439 -1.18219236 -1.3731142 -1.3986488 -1.3285032
## 6 -0.3302744 -0.3275334 -0.2247280 -0.08284935  0.1681415  0.5966643  0.7992552

#creamos el raster vacio para luego rellenarlo
result <- raster(b[[1]])

#rellenamos el raster creado con los datos del cluster del k-means
result <- setValues(result, kMeansResult$cluster)
#writeRaster(result,"./stack_fus_kmeans_1.tiff")

plot(result, col=c("red","black","yellow","green","blue","white"))

#ahora probamos la clasificacion con las imagenes normales sin modificar

b<-PES_stack_crop
b[is.na(b)]<-0
kMeansResult <- kmeans(scale(b[]), centers=6)
kMeansResult$centers

##   LC08_L1TP_200034_20230228_20230228_02_RT_B1
## 1                                   0.3236445
## 2                                   1.8669920
## 3                                  -0.8137024
## 4                                   1.0911504
## 5                                  -1.0492191
## 6                                  -0.4410001
##   LC08_L1TP_200034_20230228_20230228_02_RT_B2
## 1                                   0.3272035
## 2                                   1.8567103
## 3                                  -0.8242341
## 4                                   1.0920240
## 5                                  -1.0532609
## 6                                  -0.4167034
##   LC08_L1TP_200034_20230228_20230228_02_RT_B3
## 1                                   0.3457095
## 2                                   1.7922517
## 3                                  -0.8069011
## 4                                   1.0807066
## 5                                  -1.1028491
## 6                                  -0.2770841
##   LC08_L1TP_200034_20230228_20230228_02_RT_B4
## 1                                   0.3614099
## 2                                   1.7198005
## 3                                  -0.7065811
## 4                                   1.0632366
## 5                                  -1.1984324
## 6                                  -0.1051973
##   LC08_L1TP_200034_20230228_20230228_02_RT_B5
## 1                                   0.4536012
## 2                                   1.4626346
## 3                                  -0.4065128
## 4                                   0.9707657
## 5                                  -1.3977319
## 6                                   0.1792124
##   LC08_L1TP_200034_20230228_20230228_02_RT_B6
## 1                                   0.2807386
## 2                                   1.4214960
## 3                                  -0.2322147
## 4                                   0.8637451
## 5                                  -1.4336181
## 6                                   0.6577859
##   LC08_L1TP_200034_20230228_20230228_02_RT_B7
## 1                                   0.1151378
## 2                                   1.4539015
## 3                                  -0.1719946
## 4                                   0.7788549
## 5                                  -1.3668828
## 6                                   0.9017684

#creamos el raster vacio para luego rellenarlo
result <- raster(b[[1]])

#rellenamos el raster creado con los datos del cluster del k-means
result <- setValues(result, kMeansResult$cluster)
#writeRaster(result,"./landsat_kmeans.tiff")

plot(result, col=c("red","black","yellow","green","blue","white"))